% This function returns the moments of the data as well as the maximum negative signal for which
% a prospect theory investor would still buy. For each output, it returns the 
% point estimate as well as the confidence bounds: [estimate, low, high]

function [datamean, datastd, dataskew, privmeanstar, indiff]=bootstrap(data,delta,lambda,alpha)

confalpha = 0.05; % level of significance for confidence interval
reps = 5000;
n = size(data,2);

% point estimates
datamean(1) = sum(data)/n;
datastd(1) = (sum((data-datamean(1)).^2)/(n-1))^(1/2);
dataskew(1) = 1/n*sum((data-datamean(1)).^3)/datastd(1)^3;
% privmean is the mean private signal that would make the investor just
% indifferent between buying and selling
myfun = @(privmean)(EU_la(1/n*ones(n,1),data'+privmean,delta,alpha,lambda)-EU_la(1/n*ones(n,1),-data'-privmean,delta,alpha,lambda));
options = optimset('Display','off','TolFun',1E-10,'MaxIter',10000,'MaxFunEvals',1E6);
[privmeanstar(1),fval,exitflag,output] = fsolve(myfun,0,options);
 % indiff is the annual return the investor gets when he has a private
 % signal that makes him indifferent between buying and selling; for a
 % risk-neutral investor, it would be zero
indiff(1) = (1+datamean(1)+privmeanstar(1))^(365*5/7)-1;

% bootstrap
meanarr = zeros(1,reps);
stdarr = zeros(1,reps);
skewarr = zeros(1,reps);
privmeanstararr = zeros(1,reps);
indiffarr = zeros(1,reps);

for i=1:reps
    sample = datasample(data,n);

    meanarr(i) = sum(sample)/n;
    stdarr(i) = (sum((sample-meanarr(i)).^2)/(n-1))^(1/2);
    skewarr(i) = 1/n*sum((sample-meanarr(i)).^3)/stdarr(i)^3;
    
    myfun = @(privmean)(EU_la(1/n*ones(n,1),sample'+privmean,delta,alpha,lambda)-EU_la(1/n*ones(n,1),-sample'-privmean,delta,alpha,lambda));
    [privmeanstararr(i),fval,exitflag,output] = fsolve(myfun,0,options);
    indiffarr(i) = (1+meanarr(i)+privmeanstararr(i))^(365*5/7)-1;
end;

meanarr = sort(meanarr);
stdarr = sort(stdarr);
skewarr = sort(skewarr);
privmeanstararr = sort(privmeanstararr);
indiffarr = sort(indiffarr);

% confidence intervals

datamean(2) = meanarr(ceil(reps*confalpha/2));
datastd(2) = stdarr(ceil(reps*confalpha/2));
dataskew(2) = skewarr(ceil(reps*confalpha/2));
privmeanstar(2) = privmeanstararr(ceil(reps*confalpha/2));
indiff(2) = indiffarr(ceil(reps*confalpha/2));
datamean(3) = meanarr(floor(reps*(1-confalpha/2)));
datastd(3) = stdarr(floor(reps*(1-confalpha/2)));
dataskew(3) = skewarr(floor(reps*(1-confalpha/2)));
privmeanstar(3) = privmeanstararr(floor(reps*(1-confalpha/2)));
indiff(3) = indiffarr(floor(reps*(1-confalpha/2)));